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ABSTRACT 

(Distribution  Limitation  Statement  A) 

Four  taethods  are  surveyed  briefly  for  computing  the  probability  of  kill  (PK) 
fcr  the  offset  circle  problem  vhich  cannot  be  solved  in  closed  form.  ASA 
FORTRAN  routines  for  these  methods  are  presented.  When  ranked  according  to 
increasing  accuracy »  the  order  is  the  Wilson-Hilfsrty  (W-H)  approximation,  the 
Pcarson'-Wilscxv-HtiJferty  (P-W-H)  approximations,  the  Germond-Wcgner  (GW)  approxi- 
aatlon,  and  the  series  expansion  (S£)  approximation.  When  ranked  according  to 
increasing  machine  time  requirements^  the  order  is  C-W,  W-H,  P-W-H,  SE  with  the 
G-W  approxioaclon  requiring  30  percent  as  much  time  as  the  SE. 
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SECTION  I 
INTRODUCTION 


The  problem  at  hand  Is  to  compute  the  probability  of  kill  (PK)  of  a  point 
target  lying  on  a  plane  surface  by  a  weapon  having  a  circular  "cookie-cutter" 
lethal  area  and  delivered  with  a  circular  normal  hit  distribution  centered 
around  a  point  removed  some  distance  from  the  point  target.  Figure  1  illustrates 
the  problem.  A  point  target  is  located  at  the  origin  0.  An  aim  point  is  located 
at  A.  The  Impact  point  for  some  specific  trial  occurs  at  I.  I  is  circular 
normally  distribiited  about  A;  th?C  is,  the  coordinates  of  1  relative  to  A  are 
Independently  and  normally  distributed  v;ith  means  zero  and  standard  deviations 
ox  =>  oy  =  0.  Surrounding  the  poiiit  target  0  is  a  vulnerable  circle  C  of  radius 
R.  The  vulnerable  radius,  R,  is  a  deterministic  variable^  We  desire  to  know 
the  probability  that,  in  a  trial  yet  to  occur,  the  impact  point  I  will  fall  within 
circle  C.  Alternatively,  (but  mathematically  the  some)  we  could  have  described 
a  circle  of  lethal  radius  R  about  the  Impact  point  and  computed  its  probability 
of  enveloping  the  point  target  0. 


Figure  1.  The  Geometry  of  the  Problem 
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The  desired  probability  is 


PK  =  P  (leC;  B,  R,  o) 


dr  d0 


J 

c 

J't  exp(-r^/2<J^)dr  d0 


Equation  (1)  reduces  to 

PK  -  1/2 


[l-e-(»/of  X„  ef] 


when  B  ■  R 


PK  »  1  -exp (-R^ 720^)  when  B  »  0 


(1) 


(la) 

(lb) 


In  the  past  when  hand  methods  of  calculation  were  more  common » 
published  tables  of  PK  =  F(R/o,  B/o)  were  frequently  used.  References  1,  2,  and 
3  are  typical  of  such  tables.  Others  are  described  In  reference  4,  Alternatively » 
these  tables  (comprising  185,  15,  and  301  pages,  respectively)  can  be  transformed 
into  a  single  graph  such  as  figure  2,  which  was  first  drawn  in  the  late  1950s  by 
T.  L.  Franklin*.  Note  that  it  Is  convenient  to  make  the  plot  with  respect  to  the 
dimensionless-  coordinates  o/R  and  B/R  with  the  parameter  B/a  appearing  as  a  series 
of  radlals.  With  this  form  of  plot,  a  sink  is  seen  to  appear  at  B»»R  for  high 
values  of  R/o  (or  low  values  of  the  dispersion  relative  to  the  lethal  radius). 

To  help  explain  figure  2,  the  data  have  been  replotted  pictorially  on  figure  3. 

New  that  machine  computation  has  replaced  hand  computation  for  most  system 
analysis  studies,  it  has  become  necessary  to  push  these  tables  and  graphs  into 
the  background  and  replace  them  with  machine  routines.  This  report  presents  a 
few  useful  routines. 


*  The  data  for  figure  2  were  computed  by  PROGRAM  MOLLY  that  used  FUNCTION  PKAREA 
(inf.  series  expansion)  to  compute  PK  and  FUNCTIONS  ROSA  and  ROSE  to  interpolate 
to  desired  values  of  PK.  These  latter  functions  are  interesting  because  of 
their  accuracy  and  simplicity.  They  are  described  in  an  AFWL  technical  report 
to  be  published  as  The  Effect  of  Input  Dispersions  and  Uncertainties  on  the 
Kill  Probability  of  a  Point  Target. 
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Since  a  closed  solution*  to  this  problem  is  not  available,  we  will  have  to 
resort  to  approximations  or  a  numerical  integration  of  the  appropriate  function 
(as  was  done  in  preparing  the  referenced  tables).  Four  approximations  of  varying 
accuracy  and  speed  are  presented  in  the  following  sections. 


*  Such  a  closed  solution  exists  for  the  spherical  case.  See  Guenther  (1961) 
SIAM  Review,  Voi.  3,  No.  3,  Pages  247-250. 
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SECTION  IX 

WIXSON-HILFERTY  APPROXIMATION 


The  distributions  of  the  compovients  (referred  to  the  coordinate  axes)  of  the 
distance  from  impact  point  to  target  point  are  distributed  as  noncentral 
(Fisher  1928)  with  one-degree-of-freedom  and  noncentrality  parameter  a^. 
Following  Patnaik  (19A9)  we  can  approximate  a  weighted  sum  of  noncentral  by 
fitting  its  first  two  moments  to  an  ordinary  central  x^*  This  central  then 
may  be  transformed  to  approximate  normality  by  using  the  Wilson-Ililferty  (Vf-H) 
transformation. 

The  mean  and  variance  of  the  central  which  approximates  the  given  non- 
central  are 


m  »  1.0  +  O.S  (B/Oj)^  (2)* 

V  «  1.0  +  (B/Oj)2  (3)* 

where  »  o|  =  0.5a^;  Oj  and  being  the  component  standard  deviations  of  the 
circular  normal  delivery  pattern. 

The  equivalent  number  of  degrees  of  freedom  of  t.ie  replacement  central 
distribution  is 

n'  ^  2b^/v 

The  Wiison-Hilferty  transformation  is  then  used  to  transform  the  central  x^ 
to  approximate  normality.  This  gives  us  the  variate  t  which  is  approximately 
normally  distributed  with  zero  mean  and  unit  variance. 


*  Equations  (2)  and  (3)  are  equations  (10)  and  (11)  of  Grubbs  (Ref.  5)  as 
simplified  for  the  circular  case. 
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t  - 


(4)* 


The  technique,  as  we  have  briefly  described  it  above,  refers  to  an  offset 
circle  problem,  but  it  is  equally  applicable  (with  some  added  complications,  of 
course)  to  the  offset  sphere  problem  as  well  as  to  elliptical  and  ellipsoidal 
ncrtnal  hit  distributions.  Grubbs  (Ref.  5)  fully  describes  the  technique, 

A  suitable  translation  of  the  above  into  ASA  FORTRAN  takes  the  following 
fona; 

1  RSQ  *  ROS  *  ROS  (ROS  *  R/o)  where  now  a  *  Oj  * 

If  (RSQ)  1,  2,  3 

2  PKAREA  -  0.0 
RETURN 

3  BSQ  »  BOS  *  BOS  (BOS  «  B/o) 

W  «  1.0  +  0.5  *  BSQ 

Y  »  (1.0  +  BSQ)/(9.0  *  M  *  W) 

T  «  ((0.5  *  RSQA’*  0.333333  +  Y  -  1.0)/SQRT(Y) 

PKAREA  »  PHI(T) 

RETURN 

END 

The  normal  distribution  function  PHI(T)  is  detailed  in  the  appendix. 

To  check  out  Che  approximation  for  the  offset  circle  problem,  a  small 
program  was  run  to  prepare  a  cable  of  PK  as  a  function  of  R/o  and  B/R.  This 
information  is  presented  in  table  I  with  corresponding  values  for  the  Pearson- 
Wilsoa~Hllferty  (P-W-H)  approximation  and  for  more  accurate  approximations  as 
presented  in  the  previously  mentioned  cables.  In  each  block  Che  uppermost  value 

*  Equation  (4)  is  derived  from  equation  (20)  of  Grubbs  (Ref.  5). 
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is  the  accurate  value;  the  taiddle  value  is  given  by  the  W-8  approxltcation.  Supple- 
csenting  table  I  is  table  II  which  shows,  as  the  slddle  and  lower  values  in  each 
blocK,  the  errors  in  probability  resulting  froa  use  of  the  P-W-H  atxd  W-H  approxi- 
oations.  An  inspection  o,f  the  loventsast  values  on  table  II  (and  within  the  heavy 
lines  which  deitolt  the  probletss  for  >  10  percent)  shows  us  that  the  W-li 
spproxiaatlon  is  generaUy  rather  good  for  engineering  applicatlon.«  except  in 
the  region  near  R/o  1.0  to  1.5,  S/R  ^  1.4  were  the  errors  are  0.020  to  0.024 
in  probability*  or  -11  to  -22  percent  of  the  correct  probability.  This  is  a  bit 
ton  cauch  of  an  error  to  accept  willingly  in  a  system  analysis  study. 


*  Kathur  (Ref.  7)  indicates  that  the  cransforxation  fross  central  to 
appro'xiMte  norsality  gives  oaxixua  errors  in  probability  of  .0344.  .0x22, 

.0069  for  1,  2,  or  3  degrees-of-freedo«s.  The  equivalent  nwber  of  degrees- 
of-freedoQ  in  the  intersiediate  central  is 

n  -  2.Z/V  »  2(1  +  1/4  [B/ojV(l  +  (B/oJ^)) 

for  the  offset  circle  pvroblen.  This  equation  giv&s  n  «  2  when  55  «  0  and 
n  ~  2.8  in  the  region  of  poor  perforsante  scntlcncd  In  the  text  above.  Conse¬ 
quently,  as  Kathur  states,  the  saxisua  error  in  probability  for  the  tvn  areas 
nencioned  will  be  .012  (for  n  «  2)  and  .0D7  (for  n  *  2.8)  with  respect  to  the 
cransfoneation  iron  central  approxioace  norcsllty.  In  the  larger 

view  (the  2-stcp  transforaatlcc  foC=  ponccitral  to  approxisst^  ujnrsaiity) 
the  errors  in  probability  will  be  larger  and  as  shown  by  table  II. 
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SECTION  III 

PEARSON-WILSON-HILFERTY  APPROXIMATION 

Something  better  than  the  W-H  approximation  is  required  for  many  probability 
problems.  Pearson  (Ref.  8)  points  out  that  a  three-moment  central  chi-square 
approximation  for  the  distribution  of  a  noncentral  chi-square  may  be  quite 
accurate.  The  resulting  central  distribution  is  then  transformed  to  approxi¬ 
mate  normality  by  the  W-H  transformation.  Grubbs  (Ref.  5)  describes  this  P-W-H 
transformation*.  Applying  his  equatloiss  to  the  offset  circle  problem  we  see  that 

V'i  =  a^la^  =  1/2  (5)** 

The  third  central  moment  is 


*  Also  known  as  the  Pearson-Merrington  approximation. 

**  Equations  (5),  (6),  (7),  and  (8)  appear  in  Grubbs  (Ref.  5}  as  equations  (4), 
(21).  (5),  and  (22). 
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The  equivalent  number  of  degrees  of  freedom  of  the  central  approximation  is 


n'  -  8/6, 


(9)* 


Let  us 
tion  is 


take  4'2  =  *  l/2^”j  as  did  Grubbs.  Then  the  central  approxima- 


f 


W 


(4'2  -  +  n 


(10)* 


Next,  we  apply  the  W-H  transformation  to  approximate  normality  to  this  central 
X^  distribution  and  find  the  variate  t  which  is  approximately  normally  distri¬ 
buted  with  zero  mean  and  unit  variance. 


It  will  now  be  expedient  to  summarize  these  formulas,  remove  the  subscripts, 
and  redefine  o  as  o  » 


Let  Q  = 

Then  q  -  ^  )'  -  1/2  (f  f  -  l]  +  1 

Let  S  =  2/9n' 


K 


r 


*  Equations  (9),  (10),  and  (11)  appear  in  Grubbs  (Ref.  5)  as  equations  (24),  (25), 
and  (26). 
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u2 

Then  S  »  -j' 

36v3 

,  .  Q^/2  +  s  -  1 

and  t  *  - - 

/S 

The  quantity  t  Is  referred  to  a  table  of  cumulative  nciwl  integrals  to  find 
the  desired  probability  or  evaluated  by  a  computer  routine. 

A  suitable  translation  of  the  above  into  ASA  FORTRAN  takes  the  foliowiog  form: 

1  RSQ  »  ROS  *  ROS  (ROS  »  R/o) 

IF  (RSQ)  1,  2,  3 

2  PKAREA  «  0.0 
RETURN 

3  BSQ  -  BOS  *  BOS  (BOS  -  B/o) 

V  «  1.0  +  BSQ 

U  »  2.0  +  3.0  *  BSQ 

Q  »  U  *  (RSQ  -  BSQ  -  2.0)/(4.0  *  V  *  V)  +  1.0 
S  »  U  *  U/(36.0  *  V  *  V  *  V) 

T  «  (Q**0. 33333333  +  S  -  1.0)/SQRT(S) 

PKAREA  =  PHI(T) 

RETURN 

END 

The  normal  distribution  function  PflI(T)  is  detailed  in  the  appendix. 

As  mentioned  earlier,  tables  I  and  II  present  a  comparison  of  the  P-'W-H 
approximation  to  accurate  computations.  In  table  II,  the  center  entry  in  each 
block  gives  the  amount  by  which  the  P-W-K  approximation  errs  in  computing  PK. 

We  see  that  the  greatest  error  is  0.0143  in  probability  and  occurs  at  R/a  =>1.5 
and  B/R  *  1,0*.  Also  in  table  II,  the  upper  entry  in  each  block  gives  the 

*  Johnson  (Ref.  9)  showed  in  1959  that  Pearson's  three-moment  approximation  is 
remarkably  accurate  in  both  tails  of  the  distribution,  l.e.,  the  upper  and 
lower  5  percent  points. 
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percentage  error  in  computing  PK  by  use  of  the  P-W-H  approximation.  We  see 
that  the  greatest  percentage  error  is  -  5.1  percent  and  occurs  at  R/o  ■  1.25 
and  B/R  •»  1.2  when  PK  =  0.2;5. 

We  may  conclude  then  tliat  the  P-W-H  approximation  Is  sufficiently  accurate 
and  easy  for  most  system  studies. 
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SECTION  XV 

GESHOND-WEGNER  APPROXIMATION 


The  Gennond-Wegner  (G-W)  three-region  approximation"  generally  has  a  smaller 
percentage  error  (Table  III)  than  the  P-W-H  approximation  and  a  comparable  run 
time.  The  equations  take  the  folloirlng  form: 


PKAREA  »  1  -  ♦ 


(BOS 


-  ^R0s2  - 


(12) 


where  ^(X)  is  the  cumulative  normal  Integral  with  zero  mean  and  unit  variance. 
See  appendix. 


PKAREA  »  R0S2/A*EXP(-B0S2/A)  (13) 

where 

A  «  2  I'  R0S2/2  (13a) 

PKAREA  *  1-EXP^-ROS2/[i.416+(.397~.0159  ROS2)BOS2]^j  (lA) 


Use  equation  (12)  if 
ROS  >  3 


or 

ROS  i  1.8  and  BOS  >1.5 
Use  equation  (13)  if 

ROS  <1.8  and  BOS  >1.5 


ox 

ROS  <  1.0  and  BOS  ^  1.5 
Use  equation  (14)  otherwise. 


*  Prepared  and  submitted  by  George  Schroeter  of  the  US  Army  Aberdeen  Research  and 
Development  Center,  Aberdeen  Proving  Ground,  Maryland  21005^ 
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Equations  (12)  snd  (13)  «rt  ftofcs  GetBood  (Ref.  10),  Equation  (14)  is  froa  Wegner 
(Ref.  11),  and  ♦(x)  la  the  cunulatlve  noroal  distribution  function. 

A  suitable  translation  of  the  above  into  ASA  PCttTRAK  say  Cake  the  following 
fora: 

1  IF  (PjOS  -  l-0>  2,  2,  3 
A  «  2.0  +  ROS  *  R0S/2.i> 

2  PKAREA  -  ROS  *  50S/A  *  EX?(-BOS*EOS/A) 

GO  TO  9 

3  IF  (SOS  -  l.«)  4,  5,  6 

4  IF  (BOS  -  1.5)  7,  2,  2 

5  IF  (BOS  -  X,5)  7,  7.  3 

6  IF  (ROS  -  3.0)  5,  5,  8 

7  FKAREA  -  1.0  -  EXP  (-ROS  *  R0S/(1.41fr#- 

(0-397  -  0.0159  *  ROS  *  ROS)  *  SOS  *  SOS  **  2) 

ilO  TO  9 

8  Y  »  BOS  -  SQRT  (ROS  *  mS  -  i.O) 

PXAEEA*  1.0  -  PHI(Y) 

9  OCKJTIKDE 
RETUES 
Esa 
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SECTION  V 

A  SERIES  EXPANSION  APPROXIMATION  OF  HIGH  ACCURACY 

Although  the  approximations  previously  described  are  sufficiently  accurate 
for  system  analyses  and  are  sufficiently  easy  to  apply,  a  more  accurate  approxi¬ 
mation  to  the  value  of  the  probability  integral  may  at  times  be  desired,  even  at 
the  cost  of  more  computer  time.  A  series  expansion  approximation  is  available 
and  may  be  used  to  compute  PK  to  any  desired  level  of  accuracy. 

Vitalis  (Ref.  3)  presents  a  few  simple  equations  that  require  evaluation  to 
compute  PK.  They  are  (usinf,  the  present  notation) 


To  obtain  his  equations,  Vitalis  derived  a  conventional  Integral  (which 
cannot  be  Integrated  in  closed  form),  expanded  it  into  an  infinite  series, 
integrated  it  term  by  term,  and  thus  arrived  at  the  infinite  series  given  above. 
The  series  is  convergent  for  all  finite  and  positive  values  of  x  and  y.  The 
question  naturally  arises  as  to  how  many  terms  of  this  infinite  converging  series 
need  to  be  evaluated  in  order  to  obtain  a  reasonable  degree  of  accuracy.  Vitalis 
shows  that  the  error  involved  in  the  computations  will  be  less  than  6  x  10“^  if 
the  number  of  terms  is  T  *  5.96  *  B/o  +  3.96  for  B/o  £  3.0.  In  order  to  go  up 
to  B/o  »  16  and  to  ensure  that  the  number  of  iterations  is  T  rounded  up  to  the 
next  larger  Interger,  we  redefine  T  as 

T  -  5.0  +  6.0  *  BOS  +  0.31  *  BOS^ 
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Computational  difficulties  arise  if  the  equation  for  PK  is  evaluated  in  the 
form  in  which  Vitalis  has  given  it.  It  is  easy  to  show  that  e~y  or  e"’^  may  be 
too  small  to  be  registered  by  the  computer  or  that  y^,  x’^,  or  n!  may  be  too 
large.  For  Instance  if  B/o  =  16  then  T  =  180  and  n!  =  180!.  This  problem  is 
avoided  by  a  simple  use  of  logarithms  as  follows: 


Q(L  “  log  {"—7^ - 1  «  -X  +  m  *  log(x)  -  V  log(k) 

Le  *  m!j  k=0 

Note  that  the  bracketed  expression  of  the  form  x”/(e^n!)  represents  a 
truncated  series  expansion  of  e^  divided  by  an  infinite  series  expansion. 

An  elaboration  of  these  equations  into  a  FORTRAN  routine  follows.  This 
routine  is  effective  unless  N  is  large  as  a  result  of  B/o  being  large.  In  this 
event,  the  time  requirement  becomes  too  large.  We  handle  this  problem  by 
transferring  to  a  different  function  when  BOS  >  16. 
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FUNCTION  PKAREA  (ROS.Bl  j) 

C 

C  FUNCTION  TO  COMPUTE  THE  PROBABILITY  OF  KILL  FOR  A  WEAPON  OF  KILL 
C  RADIUS  R,  OFFSET  AIM  DISTANCE  B,  AND  CIRCULAR  NORWi  DELIVERY 
C  DISPERSION  SIGMA  (R/SIGMA=ROS,  B/SIGMA=BOS) 

C 

C  ROUTINE  BY  J.  R.  WILLIAMS,  AFWL,  KAFB,  N.  K.  87117 
C 

C  SEE  VITALIS.J.A.,  ET  AL  *TABLE  OF  CIRCULAR  NORMAL  PROBABILITIES* 

C  RPT  NO  02-949-106,  JUNE  1956,  BELL  AEROSYSTEMS  CO,  BUFFALO,  N.  Y. 

C 
C 

1  X*0.5*R0S*ROS 
XX=ALOG(X) 

IF (BOS. GT. 16.0)  GO  TO  3 
IF(BOS.LT.0.02)  GO  TO  4 
Y=0.5*BOS*BOS 
Yy=AL0G(Y) 

T=5 . 0+6 . 0*BOS+0 . 31*B0S*B0S 

N=0 

Q=-Y 

QQ=-X 

SUMTT*EXP (QQ) 

TM»EXP (Q ) * (1 . 0-SUMTT ) 

SUMT»TM 

2  N*N+1 

FNN=AL0G (FLOAT (N)) 

Q=Q+YY-FNN 
QQ=QQ+XX-FNN 
SUMTT=SUMTT+EX? (QQ) 

TM»FJCP  (Q  )  *  ( 1 . 0-SUMTT  ) 

SUMT=SUMT+TM 
IF(N.LT.T)  GO  TO  2 
PKAREA=SUMT 
RETURN 

3  PKAREA=PKAREA2 (ROS , BOS ) 

C  FUNCTION  PKAREA2  IS  THE  PEARSON-WILSON-HILFERTY  APPROXIMATION. 

RETURN 

4  PKAREA=1.0-EXF(-X) 

RETURN 

END 


L 
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SECTION  vr 

COMPUTER  TIME  REQUIREMENTS  FOR  THE  FOUR  METHODS 

A  computer  time- requirement  study  was  made  using  the  CDC-6600  computer  and  a 
program  whose  major  work  was  the  computation  of  17,340  PR  values.  Less  compiling 
time,  the  central  processor  spent  the  following  times  on  each  PK  value: 

Wilson-Hilferty  0.25  millisecond 

Pearson-Wilson-Hllferty  0.26  millisecond 

Germond-Wegner  0.12  millisecond 

Series  Expansion  0.41  millisecond,  average 

For  the  first  three  approximations,  the  time  requirements  are  essentially 
fixed.  For  the  SE  approximation,  the  time  requirement  is  a  very  strong  function 
of  BOS  (BOS  »  B/o  *  B/R  x  R/o)  and,  consequently,  increases  rapidly  as  B/R  and/or 
R/e  increase.  Thus,  the  time  requirement  may  be  either  more  or  less  than  the 
average  value  of  0.41  millisecond  Indicated  above.  If  run  on  an  IBM-7090,  the 
time  requirements  would  be  about  12  times  greater  than  Indicated  above. 

Thus,  we  see  that  the  series-expansion  method  is  the  best  and  also  the  most 
expensive.  The  G-W  method  is  seen  to  be  the  next  best  and  the  cheapest.  The 
other  two  methods  may  be  rejected  as  being  both  poorer  and  more  expensive. 
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APPENDIX 

t 

FUNCTION  PHI  (X) 

FUNCTION  TO  RETURN  THE  CUMULATIVE  NORMAL  DISTRIBUTION  FOR  A  VALUE 
OF  X  (IN  UNITS  OF  SIGMA),  WHERE  THE  MEAN  IS  ZERO. 

ROUTINE  BY  HARRY  M.  MURPHY,  JR.,  3  JULY  1969. 

AIR  FORCE  WEAPONS  LABORATORY,  KIRTLAND  AFB,  NEW  MEXICO. 

ALGORITHM  TAKEN  FROM  *HANDBOOK  OF  MATHEMATICAL  FUNCTIONS*,  (NBS 
APPLIED  MATH  SERIES  NO.  55,  1964),  SECTIONS  26.2.1  and  26.2.17, 

PAGES  931  AND  932.  (HM) 

REAL  B1,B2,B3,B4,B5,C,D 

I 

LOGICAL  PLUS 

DATA  B1,B2,B3,B4,B5,C, D/0. 319381530,-0.356563782, 1.781477937,  , 

1-1.821255978,1.33027443,0.3989422804,0.2316419/ 


AX=X 

PLUS=.rRUE. 

IF  (AX)  3,2,4 

PHI=0.5 
RETURN 

AX=-AX 
PLUS=. FALSE. 

C 

4  T=1.0/(1.0+D*AX) 

Yr^C* (T* (Bl+T* (B2+T* (B3+T* (B4+T*B5) ) ) ) ) *EXP (-0. 5*AX*AX) 
IF  (PLUS)  GO  TO  5 
C 

PHI=Y 

RETURN 

C 

5  PHI=1.0-Y 
RETURN 
END 
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